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ABSTRACT 

In one widely discussed model for the formation of nuclear star clusters (NSCs), massive globular 
clusters spiral into the center of a galaxy and merge to form the nucleus. It is now known that at 
least some NSCs coexist with supermassive black holes (SBHs); this is the case, for instance, in the 
Milky Way. In this paper, we investigate how the presence of a SMBH at the center of the Milky 
Way impacts the merger hypothesis for the formation of its NSC. Starting from a model consisting 
of a low-density nuclear stellar disk and the SMBH, we use direct A^-body simulations to follow the 
successive inspiral and merger of globular clusters. The clusters are started on circular orbits of radius 
20 pc, and their initial masses and radii are set up in such a way as to be consistent with the galactic 
tidal field at that radius. These clusters, decayed orbitally in the central region due to their large 
mass, were followed in their inspiral events; as a result, the total accumulated mass by ss 10 clusters 
is about 1.5 x 1O''M0. Each cluster is disrupted by the SMBH at a distance of roughly one parsec. 
The density profile that results after the final inspiral event is characterized by a core of roughly this 
radius, and an envelope with density that falls off p ^ . These properties are similar to those 
of the Milky Way NSC, with the exception of the core size, which in the Milky Way is somewhat 
smaller. But by continuing the evolution of the model after the final inspiral event, we find that the 
core shrinks substantially via gravitational encounters in a time (when scaled to the Milky Way) of 10 
Gyr as the stellar distribution evolves toward a Bahcall-Wolf cusp. We also show that the luminosity 
function of the Milky Way NSC is consistent with the hypothesis that 1/2 of the mass comes from old 
(~ 10 Gyr) stars, brought in by globular clusters, with the other half due to continuous star formation. 
We conclude that a model in which a large fraction of the mass of the Milky Way is due to infalling 
globular clusters is consistent with existing observational constraints. 

Subject headings: galaxies: Milky Way Galaxy- Nuclear Clusters - stellar dynamics - methods: nu- 
merical, A^-body simulations 



1. INTRODUCTION 

The centers of low-luminosity spheroids, Mb ^ — 18, 
are often marked by the presence of compact stellar nu- 
clei with half-light radii of a few parsecs and luminosi- 
ties that are ~ 20 times that of a typical g lobular clus- 
ter (iCarollo et al. lHOOl iBoker et aL]|2002t [C"6te et al. I 
120061) ^ Th e nearest such syste m is at the center of 
our galaxy (ISchodel et al. l[20Tll) . The Milky Way NSC 
is close enough tha t its radial and kinematical struc- 
ture can be resolved (jSchodel et al. |[2007t ISchodel et al. I 
l2009fl. Its total mas s is estimat ed at ~ IO^Mq 
(jLaunhardt et al. I 120021: ISchodel et al. 20 0i) and its 
half-light radius is roughly 3-5 pc Graham fc Spitler I 
|2b09; Schodcl et al. 2009, 2011.). There is a central core 
of radius ^ 0.5 pc ([Buchholz et al. II2009D . beyond which 
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the d ensity falls off roughlv as r (iGenzel et aT~ll2003l : 
Scho del et al. l[2007t lOh et al. 1 120091) . 

NSCs with properties similar to those of the Milky Way 
have now been detected in galaxies of all Hubbl e types ; 
a review of their properties is given by 'BokerJ (pOm) . 
NSCs exhibit complex star formation histories; while the 
bulk of the stars appear to always be old, the fraction 
of young stars increases toward late-type galaxies. In 
galaxies beyond the Local Group, NSCs are typically un- 
resolved, and the only structural properties that can be 
derived are half-light radii and total luminosities. The 
study of NSCs has raised considerable interest because 
of the fairly strong correlations between their masses and 
the prop erties (mass, velocity dispersion) of their hos t 
galaxies (jFerrarese et al. II2006I: IWehner fc Harris Il2006[) . 
These correlations suggest that the formation of NSCs 
and their host galaxies are linked in important ways. 
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Two models for the formation of NSCs have been 
widely discussed. In the in-situ formation model, buildup 
of molecular gas near the center o f a galaxy leads t o 
episodic star formation events fe.g.. iLoose et al. Ill982[ ). 
In this model, a NSC consi sts mostly of stars that formed 
locally (jSchinnerer et al. I [2006 . 200^ . A number of 
mechanisms have been discussed for bringing the gas 
to the center, including the magneto-rota tional instabil- 
ity in a differentially rotating gas disk (jMilosavlievic I 

f04^, tidal compression in shallow density profiles 
niscUc m fc van dc Vc n 2008)or dynam ical instabili- 
ties (jShlosman fc Begelman IH9M IBekki |[2007l ). 

Alternatively, in the merger model, globular clusters 
sink to the center of a galaxy via dyna mical friction and 
merge to form a compact stellar system ("Tremaine et al 



1975 



ICapuzzo-Dolcettanil9 93: Agarwal & Milosa vlievic 



20111) . Observations of NSCs in dwarf elliptical galax- 
ies suggest that the maj ority of such nu clei might 
have formed in this way (jLotz et al. II2OO40 . Numer- 
ical simulations have also shown that the basic prop- 
erties of NSCs are consistent with a merger origin 
JBekki et al. ''2004'; Ca puzzo-Dolcetta fc Miocchi I l200l 
jHartmann ct al. 2011). 

In addition to a NSC, the Milky Way also contains 
a massive black hole (SMBH) whose mass, M, w 4 x 
IO^Mq (jOhez et al. I l200l ICiliessen l200l . is compa- 
rable with that of the NSC. A handful of other galax- 
ies are also known to contain both a NSC and a SMBH 
(jSeth et al. Il2008t iGraham fc Spitler Il2009l ). and the ra- 
tio of SMBH to NSC mass in these galaxies is of order 
unity. The possibility of a direct link between the popu- 
lation of intermediate mass black holes that might form 
in orbitally decaying star clusters and the growth of su- 
permassive black holes at the center of galaxies h as also 
been suggested in previous papers (e.g., , Ebisuzaki et al. I 
I2OOII: iPortegies Zwart et"aLll2006[) . 

A simple argument leads to a 1 pc scale as the relevant 
one for the merger model for the formation of NSCs. The 
beginning of the disruption process of a globular cluster 
due to tidal stresses from a SMBH is expected when it 
passes within a certain distance of the galaxy center, lim- 
iting the density within that radius. Disruption occurs 
at a distance r = rdisr from the SMBH, where 
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Here p(0) is the central (core) density of the globular 
cluster, ffK its central, one-dimensional velocity disper- 
sion, and vk its core radius: the s econd relation is the 
"core- fitting formula" (jKing IIT966I ). Writing 
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for the gravitational influence radius of the SMBH, where 
cnsc is the stellar velocity dispersion in the NSC, equa- 
tion ([T]) becomes 
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Setting tk = 0.5 pc and ax = 20 km s ^, values char- 
acteristic of the most massive globular clusters, we find 



?'disr ~ 1 pc for the Milky Way. This is roughly equal to 
the radius of the core (~ 0.5 p c) that is observed in the 
distribution of late- type stars (jBuchholz et al. 1120091) . 

In this paper, we use direct A'^-body simulations to test 
the merger model for the formation of the Milky Way 
NSC. Our initial conditions consist of a SMBH and a 
diffuse stellar component that models the inner p arts of 
the nuclear stellar disk (jLaunhardt et al."1 l2002fl . The 
NSC is built up by the successive inspiral of globular 
clusters, which we inject into the system at a radius of 
20 pc. The clusters are assigned masses and radii con- 
sistent with those of globular clusters that were initially 
very massive ('^ 4 x IQ^Mq) but which were tidally lim- 
ited by the Galaxy's tidal field. As the clusters spiral in 
due to dynamical friction against the stars in the disk 
component, they eventually come within the radius of 
tidal disruption of the SMBH. We follow 12 such inspi- 
rals, resulting in the accumulation of ^ 10^ Af©. The 
NSC that results has properties that are consistent with 
the observed properties of the Milky Way NSC, including 
a p ^ density profile and a parsec-scale core. 

At the Galactic center, the relaxation time at SgrA*'s 
influence r adius is roughl y 20 — 30 Gyr, assuming Solar- 
mass stars (jMerritt Il2010l ). This is too long for a Bahcall- 
Wolf (1976) cusp to have formed over the Galaxy's 
lifetime, consistent wi t h the observed l ack of a cusp 
(iBuc hholz et"an I200I iDo et al. I I200I iBartko et aU 
2010). But a pre-existing core with radius smaller than 
the SMBH influence radius would have shrunk apprecia- 
bly o ver a time of 10 Gyr due to gravitational encoun- 
ters (jMerritt II2O1O0 . We investigate the effect of such 
evolution on our NSC model by continuing the A'^-body 
integrations after the final infall event, for a time that 
corresponds to roughly 10 Gyr after scaling to the Milky 
Way. The core radius decreases by roughly a factor of 
two in this time, bringing it to a size that is more con- 
sistent with the observed core size. The density profile 
beyond the core remains nearly unchanged. 

Since Galactic globular clusters are almost exclusively 
ancient objects with typical ages 10 — 13 Gyr (e.g., 
[Rosenberg et al. II1999I) . the merger model predicts that 
the bulk of the nuclear population is in old stars. Accord- 
ingly, using Hubble Space Telescope Near-Infrared Cam- 
era and Multiobject spectrometer (NICMOS) imaging of 
the inner 30 pc of the Galaxy, we show that the luminos- 
ity function of the Milky Way NSC is consistent with the 
hypothesis that a large fraction of its mass is in ancient 
stars. 

The paper is organized as follows. The details of our 
initial models are given in Section 2. Section 3 describes 
our simulations and results. Section 4 is devoted to the 
study of the coUisional evolution of the NSC following its 
formation. The implications of our results in the contest 
of NSCs and Galactic center dynamics are discussed in 
Section 5. Section 6 sums up. 

2. INITIAL CONDITIONS 

In the following, we describe direct N-body simulations 
that were used to study the consecutive infall and merg- 
ing of a set of 12 globular clusters each starting from a 
galactocentric distance of 20 pc. The total mass of the 
12, tidally-truncated clusters sums to ~ 1.5 x 10'' Mq, 
roughly t he observ ed mass of the Milky Way nuclear star 
cluster (|Sch6del et al. ,,2009, ). We chose 12 because this 
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number permits a particularly elegant algorithm for gen- 
erating the orbital initial conditions that does not favor 
any particular direction, as described in Section 2.2. Af- 
ter the first cluster had spiraled in to the center and 
reached a nearly steady state (as evaluated via Lagrange 
radii), we added the second cluster. This procedure was 
iterated until 12 clusters accumulated and merged in the 
inner region of the galaxy where we initially placed a 
central SMBH. Snapshots from the simulations are given 
in Figure [H This scheme has been adopted to repre- 
sent a realistic frame where the time interval between 
two consecutive globular cluster infalls and mergers to 
the center is longer than the time required for a single 
globular cluster to reach a quasi-steady state following its 
merger into the growing NSC. In any case, the results of 
the following simulations are robust and not depend (as 
far as the structure of the emerging NSC is concerned) 
much on specific initial conditions, as we discuss at the 
end of Subsection 13.11 

We now outline the details of the initial conditions 
adopted in the simulations. 

2.1. The Galactic Model 

The nuclear bulge is distinguished from the larger 
Galactic bulge (effective radius ^ 1 kpc) by its fiat disk- 
like morphology, high stellar densities, and a history of 
continuous star formation. The nuclear bulge dominates 
the inner 300 pc of the Milky Way and it appears as a, 
distinct, massive disk- like complex of stars and molec- 
ular clouds which is, on a large scale, symmetric with 
respect to the Galactic center. It consists of an r^^ nu- 
clear stellar cluster within the inner ^ 30 pc, a larger 
nuclear stellar disk and a nuclear molecular disk of same 
size (radius ~ 200 pc and scale height ~ 45 pc). The 
total stellar mass and luminosity of the nuclear bulge 
are 1.5 x 10^ Mfy^ and ^ 2.5 x 10^ Lq, respectively 
(jLaunhardt et al.Tl2002f ). The density distribution 
holds only within the NSC in the central ~ 30 pc, while, 
at larger radii, the mass distribution is dominated by the 
nuclear stellar disk which h as essentially a flat density 
profile (jSchodel et al. |[2011[ ). The initial conditions for 
the galaxy in our simulations model the nuclear stellar 
disk and they omit the central NSC. Accordingly, they 
correspond to a shallow density cusp around a SMBH, 
which is included as a massive particle, M, = 4 x IO^Mq, 
located at the origin. 

We adopted a truncated power-law model for this com- 
ponent: 



Pgx{r) = p 



n-%echf^) 



(4) 



where p = 4OOM0/pc'^ is the density at r = f = 
lOpc, and the truncation function is the same used by 

iMcMillan fc Dehnen I ()2005[ ). Since sech(a:) « 1 - ^ for 
a; ^ 1, the model is essentially a power law at r ^ Tcut, 
but it tends exponentially to zero for r 3> ^cut- We 
chose 7 = 0.5, corresponding to the shallowest power 
law consistent with an isotropic velocity distribution in 
a point-mass potential. The resulting model implies a 
mass density at 10 pc similar to what is found in the 
Galaxy outside the NSC (-- 400 Mg/pc^). We chose 
'''cut = 22 pc which gives a total mass of the (truncated) 
galactic model equal to 9.1 x 1O''M0. 



In order to generate a Monte-Carlo realization of the 
distribution function corresponding to the truncated den- 
sity profile of equation Q w e followed the method de- 
scribed in iSzell et al. I (|2005[ ). Using Eddington's for- 
mula, it can be shown that the cumulative fraction of 
stars at radius r with velocities less than v is: 
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where E = + (f>{r) and (/)(r) is the total gravitational 
potential produced by the stars and the SMBH. Once the 
positions are assigned, equation ([5]) can be numerically 
solved to distribute the particles in velocity space. 

The number of particles used to represent the galaxy 
was N — 240, 000, which implies a mass of ^ 38OM0 for 
each particle in the system. 

2.2. The Globular Cluster Model 

The globular clusters were initially placed on circular 
orbits with orbital radii tq — 20 pc. In order not to favor 
any particular direction for the inspiral, the orbital an- 
gular momenta were selected in the following way (e.g., 
iGualandris fc Merritt l[2009| ). The surface of a sphere 
can be tessellated by means of 12 regular pentagons, the 
centers of which form a regular dodecahedron inscribed 
in the sphere. The coordinates of the centers of these 
pentagons were identified with the tips of the 12 orbital 
angular momentum vectors. In this way, the inclina- 
tion and longitude of ascending node of each initial or- 
bit were determined. The choice of circular orbits was 
motivated by the well-known effe ct of orbital circulariza- 
tion due to dynamical friction (ICasertano et al. 1 119871 : 
llbata fc Lewis1ll998t [Hashimoto et al. 1120031 ). 

At a distance of 20 pc from the Galactic center, a 
globular cluster would already have been subject to tidal 
forces from the galaxy and the SMBH, and its total mass 
and radius would be less than their original values when 
the globular cluster was far from the center. We as- 
sumed that the central properties of the globular clus- 
ters were unaffected by tidal forces during the inspiral 
to 20 pc, and adopted values characteristic of massive 
clusters: central velocity dispersion aK = 35 km s~^ and 
core radius r^f = 0.5 pc. If the dimensionless central 
(King) potential is Wq = 8, the total mass works out to 
be TO ss 4 X 10^ Mq. This value of aK is roug hly two 
times the maximum value of ~ 18km listed in lHarrisI 
(|2010f )'s compilation of Galactic globular clusters proper- 
ties, while the core radius is roughly equal to the median 
value in that compilation. Our choice of such a large 
value for ax is justified by the fact that only massive 
clusters, if they are compact enough, could have arrived 
in the central regions of the Galaxy in a reasonable time 
without being destroyed by Galactic tidal forces in the 
process (Miocchi et al. 2006 and 15. 1|) . 

We then needed to generate equilibrium models for 
globular clusters with these same central properties, but 
with total masses and limiting (tidal) radii consistent 
with the known tidal forces from the Galaxy model at 
20 pc. This is not a completely straightforward exercise, 
since the gravitational force from the globular cluster act- 
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Fig. 1. — Snapshots from the A'^-body integrations, projected onto a fixed (x — y) plane at the start of each infall event. Only particles 
coming originally from the infalling clusters are displayed. The SMBH is shown as the red circle. 



ing on a star at the cluster's limiting radius, rx, depends 
both on rx and on the cluster mass mx within rx, and 
mx is a function of rx- 

We proceeded in the following way. We first assumed 
rx 3> rK- In this King-like model satisfies the 

following relation between mx and rx'- 



Gmx 
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aj^rx- 
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Given this relation, the tidal radius can then be r elated 
to th e Galactic potential and density p by (e.g.. iKingl 
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Using the galaxy mass distribution of equation (4) and 
considering the presence of the SMBH, but ignoring the 
truncation function, we find 



d<j> _ 8tt 
dr 5 
giving a limiting radius of 
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and a tidally-truncated mass from equation ([6]). Adopt- 
ing a distance ro = 20 pc we find ~ 8 pc and 



mx ~ 1.1 X 10^ M©; in other words, roughly 3/4 of 
the globular cluster mass would have been removed in 
the process of inspiralling to 20 pc. 

We then equated this mx with the mass of a new King 
model having the same core properties: 



mx = mK = p(0)r^^(Wo) 



rK^Wo)- (10) 



Here /i(Wo) is a function of the d i mens ionless central 
potential that is tabulated bv lKing 1 (|1966f) . Since all the 
quantities in equation (jlO|) are known except for Wq, we 
can solve for this variable, and find Wq = 5.8. The three 
parameters (W^Oi ''if j ) then uniquely define the King 
model that was used to generate the initial conditions of 
the globular clusters. 

As previously stated, the mass of the single particle in 
the galaxy was SSOM©. For the particles in the clusters 
we choose 2OOM0 , approximately one half of that value. 
With this choice, the total number of particles in each 
globular cluster was 5, 715 with 740 particles contained 
within the cluster core. 



Formation of the Milky Way Nuclear Star Cluster 



5 




Fig. 2. — Lagrange radii of stars from the first cluster to arrive at 
the center of the galaxy. In the upper panel the radii are computed 
with respect to the center of density of the globular cluster, while 
in the lower panel they are computed with respect to the SMBH. 
The time for each cluster to reach a steady state after its disruption 
is a few Myr. 
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Fig. 3.— Density profile of the NSC after 3, 6, 9 and 12 infall 
events. The central density grows with time. The dashed line is 
the fit to the NSC density profile obtained at the end of th e en tire 
simulation using the broken power-law model of equation lllll l . 



3. TV-BODY SIMULATIONS 

Our simulations were performed by using ^GRAPE 
([Harfst et al. 1 120071 , a direct-summation code optimized 
for computer clusters incorporating the GRA PE special- 
purpose accelerators (jMakino fc Taiii lll998D . The code 
implements a fourth-order Hermite integrator with a 
predictor-corrector scheme and hierarchical time step- 
ping. The accuracy and performance of the code are set 
by the time-step parameter rj and the smoothing length e. 
In what follows, we set r] = 0.01 and e — 0.02rK (10~^ pc 
in our case), With these choices, energy conservation was 
typically ^ 0.01% during each merging event. The simu- 
lations were carried out using the 32-node GRAPE clus- 
ter at the Rochester Institute of Technology, and also on 
computers containing Tesla C2050 graphics processing 
units at Sapienza-Universita di Roma. In the latter in- 
tegrations, 0GRAPE was used with sapporo, a CUDA 
library that emulates double- precision force calcula tions 
on single precision hardware ()Gaburov et al. |[2009[ ). 

During each infall event, we followed the evolution of 
the system until the globular cluster had reached the 
center of the galaxy and established an approximately 
steady state. This condition was verified by studying the 
time evolution of the globular cluster Lagrange radii, con- 
structed both with respect to the center of density of the 
clust er (as defined by the algorithm in Cascrtano & HuD 
I1985D . and with respect to the central SMBH. When the 
Lagrange radii had reached nearly constant values, the 
next globular cluster was introduced. Figure [5] plots the 
time evolution of Lagrange radii for the first infall event. 
The figure shows that each merging episode lasts approx- 
imately lO^yr and that the time scale for a globular clus- 
ter to reach a steady state following its disruption is in- 
deed very short, of the order of Myr. 

We evaluated the iV-dependence of our results by sim- 
ulating the first three infalls using the same orbital initial 
conditions but with ten times more particles to represent 
the clusters. Comparing the density profile of the NSC 
after the three infalls with that obtained in the origi- 
nal integrations did not reveal any significant differences 
between the two cases. 

3.1. Results: density profiles 

Figure [3] shows the density profile of the system after 
the complete merging of 3, 6, 9 and 12 clusters. We 
fitted the spatial density of the final system, within 10 pc 
arou nd the SMBH, using the broken power law model 
(e.g., [Saha1[T99l [Zhao1[l99l : 



r 

Pb 1 — 



1+1- 



{■yi-P)/c 



(U) 



where ji is the slope of the inner density profile, j3 the 
external slope and a is a parameter that defines the 
transition strength between inner and outer power laws. 
The best-fit parameters were pb — 4.1 x W^Mq/pc^, 
n = 1.5pc, 7, = 0.45, /3 = 1.90 and a = 3.73. The 
model corresponding to this set of parameters is given 
by the dashed line in Figure [3l 

Figure S] (upper panel) plots the spatial density at the 
end of the simulation over a wider radial range than in 
Figure [3l We fitted the total density as a superposition 
of two parametric models, one intended to represent the 
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Fig. 4. — Spatial (upper panel) and projected (lower panel) den- 
sity profiles at the end of the simulation. In each panel, the empty 
circles give the density profile of the A'^-body model, the solid lines 
give the best fitting model to the entire system (galaxy+NSC) and 
the dashed curves give the fit to the density profile of the galaxy, 
see text for explanation. 
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Fig. 5. — Upper panel: empty circles represent the half mass 
radius of the central NSC as a function of the total mass of the 
same syst em. The solid line represents the scaling relation given in 
equation HIGH . Lower panel: the core density of the cluster versus 
its half mass radius (empty circ les). The solid line shows the pc — r^ 
relation given by equation I I17I I. 



resented by a Sersic law: 



NSC and the other the ga laxy. For the NSC we adopted 
the modified Hubble law (|Rood et al. II1972D : 



Pci{r) = Po,ci 



r 
ro.ci 



(12) 



with best fitting parameters /9o,ci — 7.46 x 10"' Mq/ pc'^ 
and ro^ci = 1.4 pc. The galaxy remained well- fit by the 
initial "truncated" power law of equation (1), when p = 
9.91 X 10^ M©/ pc^, f = 10 pc, 7 = 0.69, and rcut = 
16.3 pc. 

The lower panel of Figure|3]shows the projected density 
of the iV-body model at the end of the simulation. We 
again fit this profile with a two-component model. For 
the NSC we used 



R 



-c 



(13) 



while the projected density profile of the galaxy was rep- 
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32.3 pc and 



The best-fit parameters were Sq ci = 2.18 x 10^ Mo/pc^ 
i?o,c/ = 1-99 pc and C = 1-03 for the NSC; Sn 

7.31 X 10''' Mq/ pc2, b = 1.68, i?0,32 

n = 1.003 for the bulge. 

Remarkably, our simulations result in a final density 
profile having nearly the same po wer-law index beyond ^ 
0.5 p c as observed (S(r) ^ r"'^: iBecklin k, Neugebauerl 
119681: iHaller et al. | [l996). In addition, the central region 
(r < rjy) of our model exhibits a shallow density profile 
(or core) near the SMBH, also in agreement with obser- 
vations (|Buchholz et al. 112001 ). The core radius in our 
model (~ 2 pc) is somewhat larger than the observed 
core, of radius ~ 0.5 pc. In 21 we show that two-body 
relaxation would cause such a core to shrink over Gyr 
time scales, as the density evolves toward, but does not 




Fig. 6. — Axial ratios of the TV-body system as a function of galactocentric radius computed after 1 (left panel), 6 (middle panel) and 12 
(right panel) infalls. Solid curves correspond to the entire model (i.e., galaxy plus the NSC), while dashed curves gives the axis ratios of 
NSC only. After the first infall the NSC is strongly triaxial in the inner regions, but appears nearly oblate at the end of the simulation. 



fully reach, a collisional steady state. 

In the upper panel of Figure [5] the half-mass radius 
(rh) of the NSC component is plotted as a function of 
the NSC mass (Md) at the end of each infall. At any 
time, the NSC mass is given by the sum of the accumu- 
lated globular cluster masses. A good fit to the data is 
obtained by 
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0.19 



pc , 



(16) 



represented by the solid line. The dependence of on 
Mci is weak, due to the fact that the size of the NSC is 
determined essentially by the fixed tidal field from the 
SMBH. 

Assuming for the growing NSC the density law of equa- 
tion ([T2|) . the core density can be defined as pc = po,ci/2^ 
and the values of pc, obtained after the end of each infall, 
can be plotted as a function of the half mass radius of 
the same system (bottom panel of Figure [5]). These data 
are well fit by 



Pc = 



1.2 X 10-^ + 1.1 exp 



1 pc 



^ . (17) 

pc^ 



shown as solid curve in the figure. 

We also carried out a separate simulation in which all 
12 clusters were placed at the same time on their initial 
orbits. Such initial conditions are hardly realistic; we 
mention the outcome briefly here since it supports the 
robustness of the results obtained from the more realistic 
initial conditions. The "contemporaneous" infall model 
produced a very similar final density profile, with a core 
of radius ^ 1 pc and a p ^ r^^ falloff at larger radii. 

3.2. Results: morphology of the NSC 

Observationally constraining the morphology and kine- 
matics of galactic nuclei is a fundamental step toward 
understanding their origin. Unfortunately, as a conse- 
quence of the strong interstellar extinction in the plane 
of the Milky Way, our knowledge of the size and mor- 
phology of the Galactic NSC are limited. Kinematic 
modeling of proper-motion data derived from the domi- 
nant population of old, mostly giant stars reveals a nearly 



spherical system of low central concentration exhibiting 
slow, approximately solid-body rotation, of amplitude 
~ 1.4 km s~Varcsec (|Trip pe et al. 2008; Schodel et aQ 
12001 . 

Aspherical NSCs are commonly observed in external 
galaxies. In a samp l e of 9 edge-on nucleated late-type 
galaxies, iSeth et alT] ()2006D reported that three of these 
galaxies (IC 5052, NGC 4206, and NGC 4244) have NSCs 
with significantly flattened isophotes and evidence for 
multiple structural components. In addition, one of these 
galaxies (NGC 4206) showed possible indication of AGN 
activity, suggesting the presence of a SMBH within the 
core of the central cluster. The NSC of the face-on galaxy 
M33, in which a SMBH is not detected (jMerritt et al. I 
120011 : [Gcbhardt et al. 2001), is also known to be elon- 
gated along an axis pa rallel to the major axis of the 
galaxy |Lauer et al. Ill998.: ,Matthews II1999I ). 

We quantified the model shape in our simula- 
tion by constructing isodensity contour s an d also 
by the momen t-of-ine rtia tensor (e.g.. iKatd 119911 : 
iPoon fc Merritt I 120041: lAntonini et al. I 12009') . as de- 
scribed in what follows: the symmetry axes are calcu- 
lated as 



Vhi/In 



T2 — \/ I22I Imax , — \/ I-^-ij Imax , 

(18) 

where la are the principal moments of the inertia tensor 
and /max = max{/ii, I22, ^33}; particles are then enclosed 
within the ellipsoid / ti'^ + y'^ / T2'^ + / t^"^ = r^. These 
previous two steps were iterated until the values of the 
axial ratios had a percentage change of less than 10""^. 
Finally, we define a > b > c letting c/a = min{ri, T2, T3} 
and b/a the intermediate value. We also define the triaxi- 
ality via the parameter T = i^a?' — 6^) / {a? — c^) . Oblate 
and prolate galaxies have T = and 1, respectively. The 
value T = 0.5 corresponds to the 'maximally triaxial- 
ity' case. 

The results are summarized in Figure |6] which displays 
the axial ratios of the NSC as a function of radius and 
at different times. The model morphology evolves from 
an initially strong triaxiality (after the first infall) into 
a more quasi-axisymmetric, oblate shape. In particular, 
the morphological structure of the final product (right 
panel) is very similar to that after the 6th infall event 
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(middle panel). This shows that the NSC is transformed 
into a nearly oblate system (T < 0.2 at r < 20 pc) af- 
ter few infalls 4), but its shape remains essentially 
unchanged from that point on. In the outer regions 
(> 20 pc), the system remained instead nearly spheri- 
cal for the entire course of the simulation. 

3.3. Results: kinematics 

Figure [7] illustrates the kinematics of the final model. 
The upper-left panel of the figure shows the one- 
dimensional velocity dispersions along, and perpendic- 
ular to, the radius vector, defined with respect to the 
SMBH. For r < 0.3 pc the system is quite isotropic, but it 
becomes tangentially anisotropic for 0.3pc < r < 20 pc. 
At larger radii the system is again roughly isotropic. The 
upper right panel of Fig. [7] shows the amplitude and ori- 
entation of the major axes of the 2d velocity ellipsoid, 
derived from the {x, y) velocities. The length of the plot- 
ted axes is proportional to the velocity dispersion. The 
tangential anisotropy is apparent. 

The lower panels of Figure [7] show the anisotropy pa- 
rameter 



2al 



(19) 



as a function of distance from the SMBH, and as a func- 
tion of time at two different distances, 10 pc and 20 pc. 
In the radial range 0.3 — 40 pc, is negative and the 
system is tangentially anisotropic. The anisotropy grows 
with the number of infall events. 

4. COLLISIONAL EVOLUTION OF THE NUCLEAR 
STAR CLUSTER 

The simulations of globular cluster inspiral described 
above took place in a short enough span of time that 
two-body relaxation effects could be ignored. T he local 
relaxation time can be defined as (jSpitzer Ill987f ) 



0.33cr3 
G^pm* In A 



(20) 



where rrii, is the stellar mass, p is the mass density and 
cr is the one-dimensional velocity dispersion. Near the 
influence radius of a SMBH, the Coulomb logarithm can 
be approximated as In A — ln(rinflcr^/2G'm*). The re- 
laxation time at the influence radius of Sgr A*, ri„fi = 
2 — 3 pc, is tr,infl ~ 20 — 30 Gyr, assuming a stellar mass 
of IMq ()Merritt |[20Tol) . This is roughly 2 x 10^ times 
the period of a circular orbit at rinfi. In our iV-body sim- 
ulations, the relaxation time is shorter (compared with 
the crossing time) by a factor of approximately 200, the 
mass of a single cluster particle in solar masses; in other 
words, it is roughly 10^ times the crossing time at Hnfl. 

In the absence of large-scale changes to the gravita- 
tional potential, an A^-body model like ours continues 
to evolve due to gravitational encounters. The evolu- 
tion that occurs should mimic the evolution that would 
take place in the real system, of much larger N ^ if the 
unit of time is taken t o be the relaxation time (e.g., 
lAarseth fc Heggie lHOOl . 

In the Milky Way, the relaxation time is short enough 
that signiflcant evolution of the stellar distribution near 
the SMBH would take place over the age of the galaxy. 
The distribution of late-type stars in the Milky Way 



NSC exhibits a near ly flat core of radius ~ 0.5 pc 
(|Buchholz et al. ir2009D . In a time of 10 Gyr, such a core 
would shrink, as the stellar density evolved toward the 
Bahcall-Woh (1976) p - r"'^!^ form inside ~ 0.2rinfl. 
Since the density profil e beyond the cor e is observed to 
have roughly this slope (|0h et al. II2009D . such evolution 
would tend to preserve the outer slope while gradually 
reducing the size of the core. A core of initial radius 1 — 2 
pc is expected to reach a size of ^ 0.5 pc, the size of the 
observed core, after ^ 10 Gyr f Merritt 1120101) . These ar- 
guments motivated us to continue the integration of our 
A^-body models after the final inspiral event. Figure [8] 
shows the density profile of the NSC at different times 
during its post-merger evolution. At the end of this inte- 
gration, i.e. after ^ 0.6ir,infl, the distribution shows an 
inner core of size ^ 0.2 pc, substantially reduced from its 
initial value of ^ 1.5 pc. The bottom panel of Figure [5] 
plots the evolution of the break radius, r&, of the best fit- 
ting broken power law profile as a function of time. The 
value of the break radius can be used as an approximate 
estimate of the model core radius. The time dependence 
of the core radius is well described by an exponential: 



r^,{t) = 1.57pcexp[-t/(0.25tr,infi)] 



(21) 



As expected, the slope of the density profile outside the 
core remains nearly unchanged during this evolution, p ~ 

^-1.8 

A core radius of ^ 0.5 pc is reached after a time of 
^ 0.25tr infi- Scaled to the Milky Way, this time would 
be 5-8'Gyr. 

Of course, in the real galaxy, it is likely that cluster 
inspiral would occur more or less continuously over the 
lifetime of the galaxy. Our separation of the evolution 
into an inspiral phase, followed by a relaxation phase, is 
artificial in this sense. Nevertheless it is reasonable to 
draw the conclusion that the size of the core resulting 
from the combined effects of cluster inspiral and relax- 
ation would be somewhat smaller than the ~ 1.5 pc that 
we found above, and therefore, closer to the observed 
core size of ^ 0.5 pc. 

Figure [9] shows the morphological evolution of the NSC 
during the relaxation phase: the radial dependence of 
the axis ratios (upper panel) and the triaxiality parame- 
ter (lower panel). There is essentially no evolution in the 
intermediate axis ratio. However, in the innermost re- 
gions of the model, the shortest axis length significantly 
increases with time. Two-body relaxation results in an 
evolution toward quasi-spherical symmetry, but at the 
end of the simulation the model has not yet reached this 
final state, still exhibiting some non-negligible triaxial- 
ity. The final model is nearly oblate with 0.3 ^ T < 0.1. 
We note that such deviations from spherical symmetry, 
although mild, might be large enough to substantially en- 
hance the rates of stellar capture and disruption by the 
SMBH with respect to the same rate computed in c olli- 
sionally resupplied loss cone theories (jMerritt II2012D . 

Figure [10] illustrates the evolution of the veloc- 
ity anisotropy profile, /3(r), during the post- merger 
phase. Relaxation tends to drive the velocity distri- 
bution toward isotropy, causing /3 to increase toward 
ze ro (this dynamical trend is sim ilar to that described 
in iCapuzzo-Dolcetta et al.1 120011) . After 0.4ti.4nfl, i.e., 
~ 10 Gyr, there remains only a small bias toward tangen- 
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tial motions, (3 ^ —0.1, r <^ 10 pc. The final anisotropy 
profile is consistent with measurements o f the Galactic 
center (jSchodel et al. |[2009HMerritt IIMol) . In the radial 
range 1" — 10", the late-type stars are observed to have a 
mean projected anisotropy of 1— {(Jt/'^r) — ~^-^'^^^i^05 
with aR and ctt the radial and tangential velocity disper- 
sions in the plane of the sky. 

5. DISCUSSION 

5.1. Cluster inspiral times 

For the model described here to be viable, globular 
clusters must spiral into the central regions of their par- 
ent galaxy in a time shorter than ^ lO^*' yr. A rough 
estimate of the time for a globular cluster of mass m at 
initial radius rg to inspiral into t he center of a galaxy as 
a result of dynamical friction is (jMerritt et al. 120041 ): 



At 



InA m \G 



(6-7)/2 



(22) 



This expression assumes circular orbits, a total mass den- 
sity of the galaxy that satisfies p{r) — pa{r/ra)~'^ , and a 
frictional force that is due entirely to stars with veloci- 
ties less than the orbital velocity of the globular cluster, 
and C oulomb logarithm InA = 6.5 ( Ant onini fc Merritt I 
120111) . The coefficient C in equation (|5.1|) is a weak func- 
tion of 7, with C = (3.9,3.6,4.2) for 7 = (1,1.5,2). 
Writing fa = ra/(lkpc), pa = Pa/ (IMqpc^^), rh = 
m / 10^ Mq and fo = ro/va, the time for a globular cluster 
to spiral in is: 



A< w 9 X lO^yr f^p^- 



?3nl/2™-lf(6-7)/2 



(23) 



Although the initial distribution of globular clusters is 
not known, a reasonable assumption is that it follows 
the distribution of the total baryonic mass predicted by 
the standard cosmological model. 

Let Atfi be the time for clusters initially within r^, the 
half mass radius of the galactic bulge, to spiral to the 
center. Within Ath, the forming nucleus has a luminosity 
comparable to the luminosities of the surviving clusters. 
An estimate of rh is re, the observed projected half-light 
(effective) radius; the stellar density at r = re is pe ~ 
14.1f~^'^^^, with re the half light radius in kpc and pe the 
stellar density in units of solar masses per cubic parsec. 
Luminosity profiles for r < re are well approximated as 
power laws with 1 < 7 < 2 in the same galaxies (Terzfc 
& Graham 2005). Thus, 



Ai w 3 X 10^° 



fl'^^rh' 



(24) 



For globular clusters of mass m ~ 10^(10®)Mo, At < 
j^Qio yj. requires re < 140(520) pc. For comparison 
bulges of SO-Sb galaxies and early type galaxies with 
-20 < Mb < -16 have 200 < re < 5000 pc , albeit 
with a large scatter (iFaber et al. lll997HFerrarese et al. I 
[200llGraham fc Worlevll2008D . Though crude, this cal- 
culation implies that a significant fraction of the globular 
clusters in faint galaxies and bulges would have spiraled 
to the center in 10^° yr. 

The dynamical evolution of the Galactic glob- 
ular cluster system has been already investigated 
by many authors usi ng different schem es and 
approximations (e.g., iMurali fc: Weinberg! 119971 : 



iTakahashi fc Portegies Zrortl [2000t iShin et all [2001) . 
None of these previous calculations were however 
designed to track the cluster orbits down to the sphere 
of influence of the central black hole. A more accurate 
estimate of orbital inspiral times for clusters in the 
Galactic bulge is obtained in what follows by adopting a 
model for the background distribution that reproduces 
more accurately the observed stellar density profile in 
the inner ~ Ikpc of the Milky Way. 

We numerically integrated the equations of motion of 
a cluster in a fixed potential including the contribution 
of dynamical friction: 



r = -V0 + adf, 



adf : 



-ATrG'^mp{r)F{< v,r) InA- 



(25a) 
(25b) 



Here p{r) is the mass density of background stars, F{< 
V, r) is the fraction of stars with local velocities less than 
that of the cluster and InA is the Coulomb logarithm. 
For the background distribution we adopted the density 
model: 



P(^) = PNSD 



PGB = 400 ^ 

pc-^ 



lOpc 



, -0.5 


r 


j exp 


70pc 



-1-6 rCxp 



(— 

VSOOpc 



(26) 



where pnsd is the density of the nuclear stellar disk and 
the second term, pgb, represents the contribution of the 
Galactic bulge which becomes significant for r > 200 pc. 
This density profi le approximates th e mass model shown 
in Figure 14 of ILaunhardt et al. I (|2002D . We chose 
m = 4 X lO^M© for the untruncated mass of the globu- 
lar cluster. In order to include the effect of tidal trun- 
cation, at any time the corresponding tidally-truncated 
mass {mx), defined in equation (jS]), was computed and 
assigned to the infalling globular cluster (i.e., m = mx) 
if mx < m. 

Figure [TT] plots the orbital evolution of globular clus- 
ters with different values of the central velocity dispersion 
and starting on circular orbits of radius 1 kpc. Globu- 
lar clusters with central velocity dispersion larger than 
about 15 kms~^ would reach the center in less than a 
Hubble time if they are initially inside the Galactic bulge. 

We stress that the calculations presented in Figure [Til 
give a conservative upper limit to the decay times for 
various reasons. Some of these are: i) the presence of gas 
in the central regions could dr amatically decr ease the 
sinking time of a stellar cluster ([Ostriker II1999D : ii) the 
background distribution is assumed to be spherical, but 
in a more realistic triaxial model the decay time could 
be greatly reduced ([Pesce et al. 1119921 ). 

In addition, the sinking times given in Figure [TT] are 
a clear overestimate of the real dynamical friction decay 
times since they were evaluated considering the Chan- 
drasekhar local approximation for tidally truncated glob- 
ular clusters on circular orbits (orbits that minimize 
the dynamical friction effect for a given orbital energy). 
More accurate orbit integrations, apt to follow the entire 
evolution of an orbit through the singularity, show that 
globular clusters on radial orbits decay in a time more 
than 30 tim es faster than on circular orbits of same en- 
ergy (jArca-Sedda fc Capuzzo-Dolcetta ,,2011.) . All this 
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Fig. 7. — Upper left: radial (cr) and tangential (at) velocity dispersions plotted as functions of the distance from the SMBH. Upper 
right: map of the principal axes of the 2D velocity ellipses derived from the {x, y) velocity components. Lower left: Anisotropy parameter 
/3. Lower right: The anisotropy parameter evaluated at 10 and 20 pc versus time. All plots were derived form the model at the end of the 
12th inspiral event. 



convince us that the decay toward the MW central re- 
gion of sufficiently massive GCs has occurred on a time 
significantly shorter than the age of the Galaxy. 

5.2. Star Formation History of the Milky Way Nuclear 

Cluster 

The NSC at the center of the Milky Way appears to 
have undergo ne continuous star format ion over the last 
10 Gyr (e.g.. iSerabvn fc Morrissi Il996f ). The evidence 
that a large fraction of the Galactic center mass was 
formed in situ at an approximately constant rate over 
the last 10 Gyr has been extensively used in the past 
to argue a gainst the merger m odel for the Mil ky Wa y 
NSC fe.g.. iMilosavlievic 1120041 : INavakshin etaiW K)^ . 

In this section, we compare the observed luminosity 
function (LF) of the stellar populations at the Galactic 
center with synthetic LFs obtained for different star for- 
mation scenarios. Our analysis suggests the possibility 
that about 1/2 of the Milky Way nuclear population con- 
sists of ^ 10 Gyr old stars brought in by infalling globular 
clusters, while the remaining mass is due to continuous 
star formation. 

Dereddened LFs of the observed populations were gen- 
erate d by Hubble Space Telescope NICMOS data taken 
from lFiger et al. I (|2004[) . The fields are all within the 



central 30 pc of the Galactic center; their exact locations 
are given in Table 1 and shown in Figure 1 of that paper. 
Data are complete at the 50% level at r7iF205W = 19.3, 
averaged over all fields. The details of the technique 
used to gene rate both model an d observed LFs are also 
described in iFiger et al. I (|2004l ) and are summarized in 
what follows. 

The LFs of the Galactic center populations are con- 
structed under the assumption that each star has the in- 
trinsic colors of a red giant and by subtracting reddening 
values for each star corresponding to its color m H — K 

or TOF160W — '71F205W- 

Geneva isochrones are used to model the LFs 
corresponding to different star formation histories. 
Th e Geneva stell ar evolu tionary mo dels ar e described 
in iShaller et al. I (I1992D . Schacrc r et al. I (:1993a, b), 
iCharbonnel et al. I (|1993l ) and iMevnet et al. I (|l994l) . We 
adopted a ma ss spectrum of stars with a Salpeter 
(|Salpeter Ill955h power law index (i.e., dN/dm = m ^■^^) 
and upper and lower mass cutoff of 120 Mq and 0.1 Mq 
respectively. The Geneva models are used to con- 
vert these masses to the absolute magnitudes in the V 
band that are subsequently transformed in the K band 
through a lookup table that relates color index to mag- 
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Fig. 8. — Post-infall evolution of the A^-body model. Up- 
per panel shows the density profile at four times: t = 
(0.12,0.24,0.36,0.48,0.6) in units of the relaxation time at the 
influence radius; line thickness increases with time. The dashed 
(blue) curve shows a fi t of the density to the broken power-law 
model of equation JTTJ at t fs! 0.36 tr.infli i.e., ~ 10 Gyr when 
scaled to the Milky Way. Lower panel plots the break (core) ra- 
dius as a function of time. The solid line is the best-fit exponential, 
equation (I21I I. 



nitude. We then sum over the histogram to produce the 
LF of each star formation event and sum the individual 
LFs to derive the LF for a given star formation scenario. 

Figure [T^] displays the results of this study for various 
star formation histories, with the model counts modified 
by the observed completeness fractions from;Figcr ct ah] 
(|1999f ). From the top to the bottom panels, the plots cor- 
respond to star formation models in which: i) the entire 
mass is build up by ancient (~ 10 Gyr old) stars, possi- 
bly brought into the Galactic center by infalling globular 
clusters; ii) the mass model is composed by ancient glob- 
ular clusters stars plus stars formed via continuous star 
formation; iii) some 1/3) o f the mass is formed durin g 
a starburst at ~ 1 Gyr (e.g., iSiouwerman et al~l 119991 ) . 
while the remainder is due to continuous star formation 
and to an ancient burst at 10 Gyr; iv) the entire mass is 
formed through continuous star formation over the last 
10 Gyr. 




r(pc) 



Fig. 9. — Evolution of the model axis ratios (upper panel) and 
triaxiality parameter (lower panel), as functions of radius, in the 
post-merger phase. Times are the same as in Figure |8] line thick- 
ness increases with time. Thickest blue lines correspond to the the 
final model. 
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Fig. 10. — Evolution of the anisotropy parameter /3 during the 
post-merger phase. Times shown are the same as in Figure [S] line 
thickness increases with time. 
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Fig. 11. — Evolution of orbital radius (upper panel) and mass 
(lower panel) for globular clusters with different central velocity 
dispersions starting on circular orbits of radius 1 kpc. These com- 
putations include the effect of tidal truncation due to the galactic 
potential (see text for explanation). Globular clusters with veloc- 
ity dispersion larger than approximately 15 kms^^ would spiral to 
the center within 10^" yr. 

The counts at faint magnitudes {Kq > 15) are con- 
trolled by ancient star formation, while the counts at the 
bright end (i^o < 8) are controlled by the extent of re- 
cent star formation activity. The brightness of the red 
clump (at Kq ^ 12) is related instead to the extent of 
intermediate age star formation activity. 

The figure shows that the ancient burst model, corre- 
sponding to a NSC composed of only ancient stars, fails 
at reproducing the observed LF. This model overesti- 
mates the counts at faint magnitudes and it does not re- 
produce the number of counts seen in the bright end. Our 
analysis rules out the possibility that the nuclear popula- 
tion consists entirely of ancient stars. On the other hand, 
star formation models in which ancient bursts are ac- 
companied by continuous star formation at other times, 
produce a LF essentially indistinguishable from that ob- 
tained when the mass is entirely due to continuous star 



formation. All these latter models fit quite well number 
counts and shape of the observed LF. We conclude that 
the observed Galactic center LF is consistent with a star 
formation history in which a large fraction of the mass 
consists of ancient ( ~ 1 Gyr) stars . A si milar result 
was obtained recently bv lPfuhl et al. I ( |2011| ). These au- 
thors constructed a complete Hertzprung Russel diagram 
of the red giant population within 1 pc from Sgr A* and 
found that about 80% of the stellar mass in these regions 
was formed more than 5 Gyr ago. 

5.3. Mass-radius relation 

In Figure [T2] the mean half-mass radius is plotted 
against total mass for nuclei (filled circles) globular clus- 
ters (open circles) and ultra compact dwarfs (UCDs, star 
symbols). We overplot the track followed by the NSC in 
our simulation during the infall events (purple-continue 
curve) and during relaxation (continue- blue line). The 
structural properties of the NSC formed in our simula- 
tions (blue-filled circle) are in good agreement with those 
of real NSCs. 

From the figure we can see that the faintest nuclei have 
roughly the same mass as a typical globular cluster. The 
size distributions for the nuclei and globular clusters also 
overlap, although the clusters in the Galaxy have half- 
mass radii of 3 pc, irrespective of mass, while the nuclei 
follow a relation of the form rh oc ^/M. Fainter than a 
few million solar masse s, the nuclei and glob ular clusters 
have comparable sizes (jHa^egan et al. II2005D . 

We now consider the merger model for nucleus forma- 
tion in the absence of a SMBH. In this case one can derive 
a simple recursive relation between the mass and radius 
of the NSC during its formation. The radius of the nu- 
cleus increases with increasing total mass, or light, as 
globular clusters merge. After the merger, its final en- 
ergy, E-f , equals the energy of the nucleus before the 
merger, Ei, plus the energy brought in by the globular 
cluster. This energy has two components: the internal 
energy or binding energy Eb, and the orbital energy just 
before the merger, Eq. From conservation of energy: 

Ef ^ E, + Eo + Eb . (27) 

Just before the merger, the orbital energy is Eg — 
aGmMi/2Ri, where Mi and Ri are the mass and ra- 
dius of the nucleus, respectively, m is the mass of the 
globular cluster, and a is a constant of order unity 
(|Hausman fc Ostriker I I1978D that depends on the ra- 
dius of the capture orbit - the radius at which the dom- 
inant influence on the trajectory of a globular cluster 
first comes from the nucleus. After the merger, the nu- 
cleus reaches a state of dynamical equilibrium quickly; 
the virial theorem implies Ef = —GM^. The equations 
above permit expressing the mass, energy, and radius of 
the nucleus recursively as 

A/,+i = (j + l)Mi, (28) 
jE,+i = ij + a)E,+jE,, (29) 

{j + if Rjl^ - ] U + a) + RY\j - 1, 2, 3, ... (30) 

where the subscript 1 denotes the initial nucleus, and, 
by assumption. Mi = m. At the time when a nucleus 
consists of few merged globular clusters, its mass and 
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Fig. 12. — Left panels display different star formation scenarios in which, from top to bottom: the entire mass is build up by ancient 
(~ 10 Gyr old) stars, that we assume to be brought into the Galactic center by infalling globular clusters; half of the mass is contributed 
by old stars and half due to continuous star formation; at the mass contributed by continuous star formation and globular clusters we 
add stars formed during a starburst episode occurred at ~ 1 Gyr; all the mass is contributed by continuous star formation. Right panels 
compare the observed LF of the Galactic center (light lines) with the LFs (heavy lines) resulting from the different star formation scenarios 
assuming Solar metallicity and canonical mass loss-rates in the Geneva models. Note that the data are much more than 50% incomplete for 
the faintest few bins. The models have not been scaled for mass, but rather have been rescaled along the vertical axis to mach the number 
counts in the K = 11.0 bin. The star formation histories corresponding to the three bottom set of panels are essentially indistinguishable 
from each other and they all reproduce quite well the observed LF. 



that of the next infalling globular cluster are compara- 
ble. In this case, equations (|28I30I) imply R cx M°-^. 
However, after many mergers, M 3> m, and the relation 
steepens to i? oc M. For a = 1.2 and 5(25)100 merg- 
ers, equations (|28l30p imply R « 2(5)10 and R cx MP, 
p = 0.5(0.6)0.7. The typic al half-mass radius o f a glob- 
ular cluster is about 3 pc ([Jordan et al. 1 120051 ) so for a 
nucleus assembled from 25 mergers, i? ~ 15 pc. This 
is in reasonable agreement with the measured sizes for 
the brighter nuclei. For a = 1.2, the expected scaling 
between rh and mass is shown by the blue dot-dashed 
curves in Figure [131 We show the predicted behavior for 
two assumptions for the mass, to, of the clusters which 
merge to form the nucleus: 10^ and 10^ Mq. At least for 
m = 10^ Mq, the agreement with the rft,-mass relation 
for real nuclei is remarkably good. 

6. SUMMARY 

We have used large-scale direct iV-body simulations to 
test the merger model for the formation of the Milky 
Way nuclear star cluster (NSC). Our initial conditions 
consisted of a massive black hole (SMBH) at the center 
of a nearly homogeneous A^-body system representing the 



nuclear stellar disk. Globular clusters were then added 
to this system, starting from circular orbits of radius 20 
pc. The clusters were tidally limited by the external field 
to have a mass of ^ 1.1 x 10^ Mq at the start. Infall was 
driven by dynamical friction, due to the stellar disk, and 
later also to the accumulated mass from the previously- 
merged clusters. The clusters were fully disrupted by the 
SMBH at a radius of approximately one parsec. After 12 
inspiral events, the accumulated mass of the NSC was 
about 1.5 X IO^Mq, comparable with the actual mass. 

The principle results of our study are summarized be- 
low. 

1 The stellar system resulting from the consecutive 
mergers has a density that falls off as ~ r~^, and a 
core of radius ^ 1 pc. These properties are similar 
to those observed in the Milky Way NSC. 

2 The morphology of the NSC evolved during the 
course of the infalls, from an early, strongly triax- 
ial shape toward a more oblate/axisymmetric shape 
near the end of the merger process. Kinematically, 
the final system is characterized by a mild tangen- 
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Fig. 13. — The measured mean half-mass radius (or effective 
radius) plotted against total mass for nuclei (filled circles) , globular 
clust ers (empty circles) an d U CDs (stars sy mbols). Data points are 
from lForbes et al. 1 1 12001 ) and lCot^ et al. | ||2^006). Blue dot-dashed 
curves show the predicted scaling in the merger model without 
SMBH for two different choices for the mass of merged clusters (see 
text for details). The red-open circle represents the initial globular 
cluster model in the A'-body simulation. Purple and continue-blue 
curves are the evolutionary track of the NSC during its formation 
and during relaxation respectively. The filled blue circle represents 
the final product of our simulation. 

tial anisotropy within the inner 30 pc and a low 
degree of rotation. 

3 In order to investigate the effect of gravitational 
encounters on the evolution of the NSC, we contin- 
ued the iV-body integrations after the final inspiral 



was complete. The core that had been created by 
the SMBH was observed to shrink by roughly a fac- 
tor of two in 10 Gyr as the stellar density evolved 
toward a Bahcall-Wolf cusp. This final core size is 
essentially identical to the size of the core observed 
at the center of the Milky Way. The density profile 
outside the core remained nearly unchanged dur- 
ing this evolution. Gravitational encounters also 
caused the NSG to evolve toward spherical symme- 
try in configuration and velocity space. 



Since Galactic globular clusters are almost exclu- 
sively ancient objects with ages ^ 10 — 13 Gyr, in 
the merger model a large fraction of the NSC mass 
is expected to be in old stars. We confronted this 
prediction with the observed luminosity function 
(LF) of the Milky Way NSC. Using stellar popu- 
lation models, we showed that the observed LF is 
consistent with a star formation history in which 
a large fraction (about 1/2) of the mass consists 
of old ( ~ 10 Gyr) stars and the remainder from 
continuous star formation. 
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